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Abstract The coronal magnetic field cannot be directly observed, but in principle 
it can be reconstructed from the comparatively well observed photospheric magnetic 
field. A popular approach uses a nonlinear force-free model. Non-magnetic forces at the 
photosphere are significant meaning the photospheric data are inconsistent with the 
force-free model, and this causes problems with the modeling (De Rosa et at, Astrophys. 
J. 696, 1780, 2009). In this paper we present a numerical implementation of the Grad- 
Rubin method for reconstructing the coronal magnetic field using a magnetostatic 
model. This model includes a pressure force and a non-zero magnetic Lorentz force. 
We demonstrate our implementation on a simple analytic test case and obtain the 
speed and numerical error scaling as a function of the grid size. 
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1. Introduction 



The solar magnetic field is observed primar ily through its effect on the polariz a tion o f 
particular Zccman sensitive spectral lines (jLandi Pegl'Innocenti and Landolfil . 120041) . 
Observations of the Fe I multiplct provide two-dimensional maps of the vector mag- 
netic field (vector magnetograms) close to the height of the photosphere. Unlike the 
photosphere, the corona lacks suitable magnetic lines and the coro nal magnetic field 
cannot be determined except under exceptional circumstances (e.g. Hvderl fl964l ) . al- 
though new m ethods based on radio and infrared observations are p resently being 
developed (e.g. IWhite and Kundulll997l [Lin. Kuhn. and Coulterl 120041) . The inability 



to observe the coronal field presents a barrier to understanding important coronal 
magnetic phenomena, including solar flares. 

In principle, the coronal magnetic field can be reconstructed from photosphcric mag- 
netograms by using a time- independent magneto-hydrodynamic model of the corona. 
The model requires solution of equations for the coronal magnetic field, B, subject 
to boundary conditions derived from vector magnctogram data. The field obtained by 
solving the model is a proxy for observational data. How accurately the model field 
reflects the true coronal field depends on the quality of the magnctogram data, and the 
accuracy of the assumptions in the model. 

No n-mag netic forces in the corona are generally negligible ( Metcalf et al. . 1199.4 
20011) and this motivates coronal magnetic field reconstructions based on a 



Gary 



force-free model, i.e. one in which only the magnetic (Lorentz) force is considered. 
In the force-f ree model the local current density is proportional to the magnetic field 
( Priest . 19841) . with the proportionality factor a in general a function of position. The 
special case where a is independent of position is called a linear force-free field. The 
linear force- free model can be solved analytically and has been extensively stu died (e.g. 
Nakaeawa and Raadull972l Barbosall978 ; Scchafcrl ll978l : lA~lissandrakisll98lf ) . However 
the linear force-free model is of limited use as a model of the coronal magnetic field 
because of specific unphysical featur es (e.g. a linear forc e- free field in an unbounded 
domain has infinite magnetic energy ( AlissandrakisL 198ll) ). 

The problems with the linear force-free model spur interest in the nonlinear force- free 
model, where a varies with position. The formal boundary value problem for force-free 
modeling can be stated in different ways depending on the choic e of boundary conditions 
on t he photosphere ( differe nt possible choices are discussed bv lGrad and Rubinl (|l958h 
and lAlv and Amarl (120071 )). Some definitions use all three components of B over the 
entire boundary, while others prescribe the normal component of the magnetic field 
B • n and the distribution of a over a single polarity of B • n, i.e. where B • n < 
or where B • n > 0. The latter has been shown to be a well-poseo3 formulation of 
the problem (in certain domains) for s ufficient l y small values of a, while the f ormer 
is an over specification of the p r oblem ( Bineau . 1972 ; Boulmezaoud and Amari . 2000t 
Kaiser. Neudert. and von Willi, l2000h . In both cases the boundary value problem is 
nonlinear and in general requires a numerical treatment. A number of different numer- 
ical solution methods have been developed for diff erent sta t emen ts of the nonlinear 
force-free bo u ndary v alue problem (see revie ws by Sakurail ( 19891 ) . or more recently 



Wieeelmannl (j2008l) ). lAlv and Amaril (j2007l) review the different methods and the 



choice of boundary conditions used in each case. 

Although appealing in its simplicity, the applicatio n of the (nonlinear) force- free 

model to magnetogram data has proved problematic (e.g. De Rosa et all^GO®. Schriiver et al^ 



A boundary value problem is said to be well-posed if it has a unique solution which depends 
continuously on the boundary conditions. 
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20081 ). The numerical solution methods are ge nerally iterativ e , and when applied to 



solar data the iterations may fail to converge ( De Rosa et al . 20091) . The problem is 



particular pronounced when large electric currents are pres ent in the data. (It sh ould 



be noted that convergence problems were not reported by De Rosa et al. ( 20091 ) for 



two of the methods which are based on a well-posed formulation of the boundary value 
problem.) 

A possible re ason f or the problems is that th e data departs significantly from a 



force- free state (|Garvl . l2001t iMetcalf et all I1995I ). This is expected to affect different 



solution methods in different ways. For methods using the well-posed formalism, two 
solutions are obtained for a single magnctogram corresponding to the two choices of 
polarity. In practice, these two solutions arc qualitatively and quantitatively different 
due to the forced nature of the boundary data. For methods which over specify the 
problem the iterative methods can never converge to a force-free state if the boundary 
data are not consistent with the force-free model. In practice, this leads to numer- 
ical solutions with residual forces and no nzero divergence. A proposed solut i on to 
this problem is 'preprocessing' of the d ata ( Wiegelmann. Inhester. and Sakurai . 2006L 



Fuhrmann. Seehafer. and Valoril 120071 ). With this technique the data are modified to 



be consistent with small net forces and torques in the overlying volume while main- 
taining minimal departure from the original data. Howe ver, the prepro c essed data 



are not necessarily consistent with the force- free model ( De Rosa et al . 20091) and 



mixed re sults have been obtai ned in the modelling: in some cases the problems were 



reduced ( Schriiver et all 120081 ) and in other cases there was no significant improvement 



(|De Rosa et al.l . l2009h . 

Another approach is to develop static models of the coronal magnetic field which 
incorporate non-magnetic forces, i.e. magneto-hydrostatic models, which incorporate 
gravity and gas pressure forces. The special case of a model with only magnetic and 
pressure forces may be called a magnetostatic model. This model requires the pre- 
scription of pressure as a boundary condition. Presently, the photosphcric pressure 
distribution is difficult to obtain observationally compared to the magnetic field. This 
is a limitation of magnetostatic (and magneto-hydrodynamic) models in application to 
solar data. 

A number of methods for solving the magneto-hydrostatic and magnetostatic equa- 
tions have been developed for modeling the coronal magnetic field. The optimization 
procedure developed for the nonlinear force-free model has b een extended to magneto- 
hydrostatic and magnetostatic models in Cartesian geometry ( Wiegelmann and Inhesterl .[ 



2003 : Wiegelmann and NeukirchL 20061) . and sp herical geometry ( Wiegelmann et al 
2007 ). Amari. Boulbe. and Boulmezaoud (I2009T) presented a finite element imple: 



men- 



tation of the Grad- Rubin method ( Grad and Rubin . 1958) for solving a magnetostatic 
model which is applicable in arbitrary geometry. In addition semi-analytic approaches 
have been deve loped but these provide only restricted, and not general, solutions (e.g. 



Rudenkoll2001l ). 



In this paper we present a numerical code to solve a general magnetostatic model 
of the coronal magnetic field of an active r egion . Our code is an implementation of 



the Grad-Rubin method (jGrad and Rubin 
the half-space 



19581 ) to the magnetostatic equations in 
The Grad-Rubin method has previously been used to solve 



z > 

the nonlinear force- free equations in the half-space for coronal reconstructions (e.g 
Wheatland! l2006t lAmari. Boulmezaoud. and Mikic 1999 ), and to so lve the magneto 



static equ ations in toroidal geometry dReiman and Greenside . 1986 ) and in arbitrary 
geometry ( Amari. Boulbe. and Bouhnezaoud 20091 ). The Grad-Rubin method has the 
advantage of, in principle, solving the magnetostatic equations to the limited imposed 
by numerical accuracy (e.g. truncation error due to the finite grid size), assuming the 
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Grad-Rubin iteration procedure converges. We apply the code to a simple analytic test 
case as a demonstration, and to investigate the speed and numerical error scaling of 
the implementation. The resulting code and method does not yet provide a practical 
tool for coronal field modeling from solar data, because of the neglect of the gravity 
force, but it is a significant step in this direction. 

This paper is structured as follows. Section [5] presents the magnetostatic equations, 
and the boundary value problem to be solved. Section [3] is a brief summary of the 
Grad-Rubin method, and describes our specific implementation of the method in code. 
Section|4]presents the analytic test case used. Section[5]prcsents the results, and Section 
[6] contains discussion of the results and the conclusion. 
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2. Magnetostatic Equations and Boundary Value Problem 

2.1. The Model 

The magneto-hydrostatic equations with the gravity force neglected are 

V x B = imjJ, (1) 

V-B = 0, (2) 

and 

J x B - Vp = (3) 



(jPriestl . 119841) . Here p is the gas pressure, J is the electric current density, and B is the 



magnetic field vector. 

2.2. The Boundary Value Problem 

For a localized active region the curvature of the photosphere is small, so we solve 
Equations (H}-© m the half-space z > 0, with the z = plane representing the 
photosphere. The appropriate boundary conditions prescribed on the z = plane are 
B z , together with p and J z prescribed over a single polarity of the magnetic field, 
i.e. p and J z are presc ribed only at points with B z > or at points where B z < 



(jGrad and Rub in. 1958). We denote the boundary values of the magnetic field, pressure 
and current density by -B bs , Pobs and Jobs respectively. These boundary conditions are 
believed to be the correct physical boundary conditions for the magnetostatic model. 
Equations (|T|)- (J3J) arc solved numerically in the finite volume 

n b = {(x,y,z)\0< x<L Xl 0<y<L y ,0<z< L z }. (4) 

In addition to the z — plane, Qb has five additional plane boundaries on which 
boundary conditions are required. Magnetogram data only provide boundary conditions 
on the z — plane and reasonable assumptions must be made for the remaining five. 
This problem is faced by all reconstruction codes regardless of the particular model 
or method used, and models including more physics typically require more boundary 
conditions at each boundary. 

We choose boundary conditions on the magnetic field such that all field lines are 
connected to the lower boundary at two points (i.e. there are no open field lines). The 
need for this is discussed in Section l5.2.1l We achieve this in practice by imposing either 
(i) closed boundary conditions on the top and side boundaries 

Bn = 0, (5) 

where n denotes the unit normal vector to each boundary, or (ii) a closed top boundary 
condition 

B • z = 0, (6) 

together with periodic boundary conditions on the side boundaries. For the periodic 
case any field line which leaves the computational volume by a side boundary re-enters 
on the opposite side, and therefore eventually connects to the lower boundary. These 
boundary conditions specify that magnetic flux only enters and leaves through the 
lower boundary, which requires that the lower boundary is flux balanced. 
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3. Numerical Implementation of the Grad-Rubin Method 

In this section we outline our implementation of the Grad-Rub in method in code . The 
approach is similar to that for the force- free code described in I Wheatland! (|2007l ) . 

3.1. The Grad-Rubin Method 

The Grad-Rubin method (also called the curre nt-field iteration method ) is an iterative 
scheme for solving the magnetostatic equations ( Grad and Rubin1 . ll958l ). In this method 



the nonlinear equations (Equations ([Tj ) -([3" ]) ) are replaced with a set of linear equations 
which are solved at each iteration. 

Here we briefly outline a sing le Grad-Rubin iteration (a more detailed description is 



given bv lGrad and Rubinl (|1958|) ). We denote a quantity after k Grad-Rubin iterations 
using a superscript, so for example B( fe ) denotes the magnetic field in our computational 
volume after k iterations starting from an initial magnetic field B^ ). In practice the 
iteration is initiated with a potential field in the volume calculated from £> b s , which 
we denote as B'°) = Bo- In the following we assume B^* 1 is known from a previous 
iteration, or is the initial potential field. A single iteration consists of the following 
steps. 

1. Calculate a new pressure p( k+i ) in the volume by solving 

B « . Vp {k+1) = 0, (7) 



with boundary conditions 



P (k+1) \ z =o=Po bs (8) 



prescribed over one polarity of -Bobs- 

2. Calculate the component of the current density perpendicular to the magnetic field 
in the volume using 

J { 1 +1) = Vp (fe+1 > xB( fe )/|B( fe )| 2 . (9) 

3. Calculate the component J|| fc of the current density parallel to the magnetic field 
in the volume. The parallel component can be written as 

j[ fc+1) = ir^BW/tt,, (10) 

where a^ k+1 ^ is a scala r function of position. T he parameter o-( fe+1 ) is calculated in 
the volume by solving ( Grad and Rubin . 19581) 



V(T (fc+i). B W = _ MoV .j( L fe+i) ; 
with boundary conditions 

°U=o = °itt 1} (12) 

where 

r(fe+l) 



(fc+1) _ Mo(J± ' g ~ ^obs) 

obs — r> 

ts, 



Sobs 



(13) 

z=0 

The boundary conditions cr^^ are prescribed over a single polarity of -Bobs and 
are calculated using Equation (fT3"j) at each iteration. The form of Equation (Tl3"|) is 
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such that j( fe+1 ) • z = J b s on the lower boundary (over the polarity of B a b s chosen 
for the boundary conditions). 

Equations ((9|)- ([T0|) define the total current density in the volume: 

j( *+i) = j(fc+D + tr (*+i) B (k)/ W) . (14) 

4. Calculate the new magnetic field in the volume by solving Ampere's law 

VxB( k+1 )^„J (W) , (15) 
where the boundary conditions on B' fc+1 ) are those in Section [21 



3.2. Overview of the Code 



The code solves the magnctostatic equations in the finite Cartesian domain fib using 
the Grad-Rubin method. The numerical grid is uniformly spaced with N grid points 
along each dimension, so there are N 3 grid points in total, and in the following we write 
L = L x = L y = L z for simplicity. The grid spacing is L/(N — 1), and the coordinates 
of the grid points are (xi, yj 7 Zk) = k)L/ (N — 1), where < i, j, k < N — 1 . 

In the following wc describe the implementation of the Grad-Rubin steps identified 
in Section 13.11 in detail, identifying some of the numerical methods used. The code 
i s an implementat ion in FORTRAN90 using double precision floating point numbers 
19921 ). The code is parallelized for shared memory multiprocessors using 
1990h . 



(Press et al 



OpenMP ([Chandra et al 



3.2.1. Step 1: Update of Pressure in the Volume 



Equation ([7]) is solved by a field line tracing (or characteristic) method. The same 
method is used in so me existing force-free Grad-Rubin method implementations fo r 
solving B-Va = fe.o. lAmari. Boulmezaoud. and Mikidll999t rWheatlandl200fi[ 120071) . 
where a is the force-free parameter defined by V x B = aB. Equation has also 



been solved using a finite element method (|Amari. Boulmezaoud. and Ahl 120061 ). In 
the force-free case the tracing is used to update a in the volume, and here it is used to 
update the pressure. 

The procedure is as follows. For each grid point (xi,yj,Zk) the field line threading 
(xi,yj,Zk) is traced to the point (xo,yo) where it crosses the lower boundary with the 
appropriate polarity for the boundary conditions. The pressure at the grid point is then 
assigned to be equal to the boundary value: 



p(xi,Vj,Zk) = Pobs(x ,y ). 



(16) 



This procedure solves Equation ([7]) because the pressure is constant along a magnetic 
field line. The tracing is performed in either the forward or backward direction along 
the field line, with the direction chosen such that the point (xo,yo) has the appro- 
priate polarity for i? b s - Fourth order Rungc-Kutta is used for the numerical tracing 



([Press et all 119921) . Because the path of the field line is not confi ned to grid points , 



trilinear interpolation is used to estimate B between grid points ([Press et a?.ITl992|) . 
Similarly, bilinear interpolation is used to compute p bs at the boundary point (xo, yo), 
which also may not coincide with a grid point. 

Open field lines (see Section ^. 2[) would present a problem for this method. Since p bs 
is only prescribed over a single polarity, it would be impossible to assign pressure to 
points threaded by field lines which connect to the lower boundary at only one polarity, 
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opposite to that for which p bs is prescribed. We avoid this problem by preventing open 
field lines through the choice of boundary conditions on B explained in Section 12.21 
This eliminates the problem but introduces artificial boundary conditions on the top 
and side boundaries. However, the region of interest can be isolated from the effects of 
the boundaries by using a large domain. 

3.2.2. Steps 2 and 3: Update of the Current Density in the Volume 

Given the updated values for the pressure in the volume, Equation (jH]) may 

be directly evaluated at every gridpoint to give the perpendicular current Jj_ in the 
volume. The derivatives in the gradient on t he right hand side are evaluated numerically 



using a centered difference approximation (IPress et all 119921) 



A particular problem is encountered with the evaluation of the perpendicular current 
J ± in the volume at each iteration (step 2 in the enumeration of the procedure in Section 
13. ip . for the test cases considered here. The perpendicular current is calculated using 
Equation ©: 

J^ +1) =Vp {k+1) xB^/|B (fc) | 2 . (17) 

For the analytic solutions we use in Section [5] there are locations in the volume where 
B = and is finite. The perpendicular current is not correctly evaluated numerically 
at these points. To prevent this we choose grid sizes for our problems such that the 
points with B = fall between grid points. This removes the problem, which is due to 
the artificial nature of the test cases. 

The value of cr( fc+1 ) at each gridpoint is then obtained by solving Equation (ITT1) . 
This equation may be integrated along a field line to give the formal solution 

^ k+l \x uyj ,z k )=<7 ohs {xo,yo)-l f V J^ +1) (x( S ))/|B( fc >(x( s ))M s , (18) 
where 



_ f +1 if Jobs is prescribed over B z > ... Q , 

^ \ —1 if J bs is prescribed over B z < 0, 

where x(s) is the path of the field line (the parameter s is the arc length along the 
field line), and where x(so) is the point (xo, yo, 0) at which the boundary conditions on 
CT (fc+i) are imposed. Trilinear interpolation is used to assign values in the argument of 
the integral along a field line, and the integral is evaluated using the trapezoidal rule 



( Press et all 1992) 



The field line tracing needed for steps one and two (updating pressure and a) is 
the computationally slow part of this implementation of the Grad- Rubin method. The 
number of operatio ns for the field line tracing scales as N 4 for a grid with N 3 points 
(|Wheatlandl . l2006h . In the following we will write '~ iV 4 ' to denote such a scaling. 



The code parallelizes the process using OpenMP, with the workload divided such that 
different code threads trace different field lines. 

It is important to understand the accuracy of the numerical solutions for p( k > and 
<j( fc ) _ p r the field line tracing solution we can infer an approximate scaling fo r the nu 



merical error as follows. Trilinear interpolation has truncation error ~ 1/N 2 (jZikanovl . 



2010), and on average a field line requires ~ N Rungc-Kutta steps to reach the lower 



boundary. Therefore the total numerical error introduced by tracing a field line to 
the lower boundary has scaling ~ N x 1/N 2 = 1/N). We expect this error to be the 
dominant error in the calculation, and so the error scaling for the whole computation 
is ~ 1/N. This error scaling has been confirmed for the force-free case ( Wheatland! . 
20061) . 
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3.2.3. Step 4-' Update of the Magnetic Field in the Volume 

The magnetic field may be expressed as B^- 1 = Bo + Be , where Bo is the potential 
field satisfying V x Bo = together with the boundary condition 

B • z| 2=0 = Sobs, (20) 

and where B c ^ is a current carrying field satisfying V x B c *^ = hqJ^ together with 
the homogeneous boundary condition 

B^ • z| z=0 = 0. (21) 

(k) 

The fields Be and Bo have the same boundary conditions on the top and side bound- 
aries (either closed boundary conditions on all other boundaries, or closed top boundary 
conditions and periodic side boundary conditions, as discussed in Section |2T2|) . 
For the potential field a scalar potential <j> can be introduced defined by 

V0 = -B , (22) 

and the problem reduces to solving Laplace's equation (j Jackson , Il998t) : 



V 2 = 0. (23) 



Lapla ce's equation has well-known Fourier solutions in Cartesian geometry ( Morse and Feshbachl .[ 



1953|), and we use a two-dimensional Fourier series solution with the expansion per- 
formed in the x and y directions. For periodic boundary conditions the components of 
our potential field ar^| 



B 0x (x,y,z) = J2Y, CmnlkmCOsh W z - L V et{Xkm+ykn) > ( 24 ) 

m—0 n—0 
oo oo 

B Qy (x,y,z) = ^^w^cosh(/c[^-L])e^ fe - + ^), (25) 



and 



rn— n—0 



oo oo 



B 0z (x,y,z) = J2 ^c mn fcsinh(fc[z-L])e I ^ + ^"). (26) 



m— n—0 

„2 _ ;„2 i 1,2 



where k m = 2-Km/L, k n = 2nn/L and k = kf n + k^. The Fourier series coefficients are 
derived from the boundary conditions: 



f\J o J a dxdyB ohs ^,y)e- t(xkm+ykn) /smh(kL). (27) 



As explained in Section 12. 2[ we also use a solution with closed boundaries, and the 
expression for the components of this field are similar (see Appendix A). Equations 
(|M]) - ([23)l can be computed using Fast Fourier transforms, in which case ~ N 3 \og(N) 



2 We note that this solution is a special case of the linear force-free solution due to Barbosa (1978), 
and is obtained by setting a = in that solution. 
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operations are req uired to evaluate the potential field for a grid of iV 3 grid points 



([Press et all 119921) . This makes step 4 relatively fast, computationally. 

For the non-potential component B c we use a three-dimensional Fourier series so- 
lution to the vector Poisson equation, wo rking with a vector potential A c such that 



V x A c (jMorse and FeshbachL 119531 ). The components of the field arc 



oo oo oo 



B C x(%,y-> z ) — ^] ^] ^] kp a rnnp 



and 



m— n— p— 

oo oo oo 

Bcy(x, y, z) = [^p fl mlp ~~ k 

m— n— p—0 



OO OO OG 



ia (3) 

m ,u itmp 



cos^pz)^*™***"")/* 2 , (28) 

cos^zK^+^V^, (29) 

(30) 



B CZ 0W) = E EE 8 [ k ^n P - k n^n P ] Sm(k pZ )e^ +k ^ / k" . (31) 
m— n— p—0 

Here fc m = 2irm/L, k n — 2nn/L, and fc p = irp/L. The coefficients amnp, with z = 1,2,3 
are given by: 



(1) 2^0 / L / £ / L 
mnp " £ 3 7o io 7o 



2/io 



and 



"mnp jj^ 



(3) = ^MO 

mnp 



o Jo Jo 



J x (x, y, z ) c -^ k ^+k n y) s [ n (k pZ )dxdydz, (32) 
J tf (s, y, z )e- li - k ™ x+k ^ sm(k p z)dxdydz, (33) 



JO JO 



J z (x, y, z)e i ( fc "» a,+fc »i') cos(k p z)dxdydz. 



(34) 



The coefficients given by Equations (j3"2")) - (|3"4")) . and the solution given by Equations i 
([3T]) . arc computed in ~ TV 3 log(T V) operations usin g a combination of fast Fourier, fast 
sine, and fast cosine transforms (|Poularikasill99d) . The corresponding solution for B c 
with closed side boundaries is given in Appendix B, which also provides more detail on 
how these solutions are derived. 
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4. Analytic Test Case 



To test the code w e use a simple analytic so l ution to the magnetostatic equations 
( Wiegelmann . 19981 Wiegelmann and Inhester , 2003 ) which is a generalization of a 
force- free field in ISturrockl (|1994l ). and which may be derived using the generating 
function method. The solution describes a sheared magnetic arcade with translational 
symmetry in the y direction and periodicity in the x direction. 

For the problem at hand we modify the sheared arcade solution by imposing a closed 
top boundary condition (to match the boundary condition required by our Grad- Rubin 
method) : 



B z (x,y,L) = 0. 
The components of the resulting magnetostatic field are 

B x {x,z) = ifj lsm(kx) cosh\l(L — z)], 
B y (x,z) = tpoXy/l — ao sin(fca;) sinh[Z(X — z)], 



and 

B z [x, z) = ipok cos(fcx) sinh[/(i — z)], 
and the pressure is given by 

p(x, z) = ipQ ( — — J sin 2 (fca;) sinh 2 ^(z — L)} 



(35) 

(36) 
(37) 

(38) 
(39) 



where k, A, ipo and ao are free parameters subject to < ao < 1 and I = \/k 2 — A 2 . 

The parameter A determines the currents in the volume, with A = giving a potential 
field. The parameter ao sets the pressure in the volume, and the special case ao = is 
a linear force-free field with force-free parameter a = A. The parameter ao may be used 
to enforce closed side boundary conditions in y with the choice ao = 1. The parameter 
k determines the period of the solution in the x direction, and the constant which 
determines the magnitude of the field, is chosen to be 



1 



sinh(ZL)fc ' 



(40) 



to specify max(|_B b s |) = 1 on the lower boundary. 

The boundary conditions in the z = plane with the choice of ipo given by Equation 
flUD are 



B obs = cos(fcx), (41) 

Pobs = o" ^ ( fcx )' ( 42 ) 

and 

Jobs = Mo-Vl - ao cos (kx). (43) 

A schematic diagram of the field lines of the solution is shown in Figure [TJ for the 
choices k — 2ir/L, ao — 0.5, and A = tt/ (2L). The view in the left panel of the Figure is 
along the y axis. This perspective shows the arcade-like field line structure. The right 
panel of Figure [1] shows a top-down view of the particular field lines shown as dashed 
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Figure 1 . Field lines for the analytic magnctostatic field test case (see Section |4j with parameters 
k = 2n/L, ao = 0.5, and A = tt/(2L). In the left panel the point of view is along the y axis. The right 
panel shows a top-down view of the field lines indicated by the dashed curves in the left panel. These 
field lines are sheared with respect to the x axis, at an angle of fa 10° . 



curveswh in the left panel. The top of these field lines is at z = 3L/A. This perspective 
shows that the arcade is sheared. For the given solution the shear angle is sa 10° for 
the field lines with the height shown. 

5. Results 

In this section we apply the code to the test case in Section @] for two different choices of 
parameters. The first choice is for a calculation with periodic side boundary conditions, 
and the second is for a calculation with closed side boundary conditions. For both cases 
the tests are performed several times on grids of varying size. The convergence of the 
Grad-Rubin iteration is demonstrated, and the speed of the code as a function of 
problem size (introduced in Section 13.2.2)) is confirmed. We also examine the accuracy 
of the numerical solution by comparison with the analytic solution, and investigate the 
accuracy as a function of grid size, for comparison with the estimate of the scaling of 
the accuracy given in Section f3. 2.21 

5.1. Test Case with Periodic Side Boundaries 

For the first test we use the parameters k = 2ir(l — 1/N)/L, A = n/ (2L) and ao = 0.5. 
The parameter k is chosen to vary with iV so that the periodicity of the solution 
matches the periodicity of the discrete Fourier transform |Press et all \l992h . A side 
effect of this is that the test case differs between grid sizes. In particular the boundary 
conditions, which depend on k through Equations (|4lj) - P3")) . are different. However, for 
large N the difference is small. We apply the code to the test case with three different 
grid sizes (N = 51, 101, 151), and in each case we perform 30 Grad-Rubin iterations 
starting from a potential field. 
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Figure 2. The mean absolute change in the field AS avg as a function of iteration number for the 
application of the Grad-Rubin method to the test cases. The left panel shows AB avg for the three 
tests with periodic side boundary conditions, and the right panel shows A_B avg for the three tests with 
closed boundary conditions. The different symbols represent AB avg for different grid sizes N. In both 
panels the N = 65 case is shown with squares, the JV = 101 case is shown with diamonds, and the 
N = 151 case is shown with plus signs, and the scale on the y axis is logarithmic. 



The convergence of the Grad-Rubin iteration is measured by the absolute change in 
the field at each iteration, i.e. 

A^gHIB^-B^- 1 )!) (44) 

where () denotes the average over all points in the computational volume. The left 
panel of Figure [5] shows A£? avg over 30 Grad-Rubin iterations for the three different 
grid sizes. The squares show the case with TV = 65, the diamonds show the case with 
N = 101, and the plus signs show the case with N = 151. In all three cases the Grad- 
Rubin iteration converges: Ai3 avg decreases exponentially for rj 15 iterations before 
becoming approximately constant. The rate of convergence does not appear to depend 
strongly on grid size. 

The numerical solution after 30 iterations is compared to the analytic solution. 
This is done qualitatively by comparing the field lines of the analytic solution, the 
initial potential field, and the numerical solution. The left panel of Figure [3] shows the 
field lines of the analytic solution (red field lines), and the field lines of the potential 
field (blue field lines) viewed looking down on the computational domain. There is 
a significant difference between the two sets of field lines. The right panel of Figure 
[3] shows the field lines of the analytic solution (red field lines), and the field lines of 
the numerical solution after 30 Grad-Rubin iterations (blue field lines) from the same 
viewpoint. The two sets of field lines closely coincide, confirming that the Grad-Rubin 
iteration has converged to the analytic solution. 

In addition to this qualitative comparisons, we compare the analytic and numerical 
solutions quantitatively. Several metrics have been devel o ped f or this purpose in the 



context of nonlinear force- free modeling ( Schriiver et ai . 20061) . The various metrics 



show similar results, so for brevity we present only the mean vector error 
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Figure 3. Comparisons between the field lines for the analytic solution, the Grad-Rubin solution and 
the initial potential field for the test case calculation in Section |5. II The left panel shows the analytic 
solution (red field lines) and the potential field (blue field lines). The right panel shows the analytic 
solution (red field lines) and the numerical solution after 30 Grad-Rubin iterations (blue field lines). 
The solutions are viewed looking down from the top of the computational domain. The z = plane is 
shaded to show -B D b s (regions with B a ^ s < arc dark and those with B ^ s > are light). The solution 
after 30 Grad-Rubin iterations closely matches the analytic solution. 



*- = (V>- ™ 

and the metric 

£ <= s = 1 -(™)' (46) 

which is based on the Cauchy-Schwartz inequality. In both these definitions the ana- 
lytic solution is B, the numerical solution is b, and the average is over all grid points 
in the domain. For two exactly matching fields E m — and Ecs = 0. The mean vector 
error is a measure of the difference between the magnitude and direction of B and b 
at each grid point, while Eqs is only sensitive to the differences in the direction of the 
fields. 

We also compute 

£div = <|V-b|>, (47) 

which is a measure of the divergence of the numerical solution b. A second-order 
finite difference ap proximation is used to compute the divergence at each grid point 



(jPress et all 119921 ). In principle -Ediv should be zero, but this is not achieved in prac- 



tice due to the finite numerical accuracy of the solution, and the truncation error 
introduced by the numerical approxi mation to differen tiation. The truncation error in 



the derivative has a scaling ~ 1/N 2 (jPress et all 119921) 



3 These metrics are chosen because they show the largest discrepancy between the numerical and the 
analytic fields. 
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Figure 4. Error metrics Em (see Equation 1451 ) and Ecs ( see Equation 1461 1 ) as a function of grid 
size N for the initial potential field (plus signs and crosses), and for the numerical solution (diamonds 
and squares). The left panel shows the results for the test cases with periodic side boundary conditions 
after 30 Grad-Rubin iterations, and the right panel shows the results for the test cases with closed 
side boundary conditions after 50 Grad-Rubin iterations. The Grad-Rubin solutions show power-law 
scaling Em ~ N 1 and Eqs ~ N"< , and the dashed lines show power-law fits to the data. The power-law 
indices 7 for each fit are summarized in Table [T] 



1000 



The left panel of Figure 0] shows E m and Ecs for the initial potential field (plus signs 
and crosses), and for the numerical solution after 30 Grad-Rubin iterations (diamonds 
and squares), as functions of grid size N. The error associated with the potential 
field is independent of N. For the Grad-Rubin solution both E m and Ecs decrease 
approximately as power laws, i.e. show scalings ~ TV 7 , where we find (based on least 
squares fits) 7 = —1.3 for E m and 7 = —2.3 for Ecs- The power-law fits are shown in 
Figure E] by the dashed lines. The scaling for E m is close to the expected ~ l/N scaling 
discussed in Section f3. 2. 21 while the scaling for Ecs is closer to ~ 1/N 2 . 

The left panel of Figure [5] shows E^iv as a function of N for the numerical solution 
after 30 Grad-Rubin iterations (diamonds) and the analytic solution (plus signs). The 
the two data sets overlap indicating that the divergence of the numerical solution is 
close to or smaller than the truncation error in the numerical derivative. In common 
with E m and Ecs t the metric -Ediv for the Grad-Rubin solution has power- law scaling 
~ iV 7 . We estimate 7 = —2.0 based on a least squares fit to the data for the Grad-Rubin 
solution (the power- law fit is shown as a dashed line in the figure). 

We also measure the run time of the code for different grid sizes using the CPU clock 
(the tests are performed on an eight core CPU). Figure [6] shows the execution time for 
30 Grad-Rubin iterations as a function of N. The results appear to follow a power law 
~ A^ 7 , and we estimate 7 = 3.8 based on a least squares fit to the three data points 
(the fit is shown by the dashed line) . This scaling is close to the ~ iV 4 scaling expected 
for the field line tracing discussed in Section f3.2.2l This scaling implies that the field 
line tracing is the computationally slowest step in the calculation. 
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Figure 5. The absolute value of the divergence of the numerical solution and initial potential fields, 
averaged over the computational domain -Bdiv ( see Equation I47H ) as a function of grid size N. The left 
panel shows E^i v for the analytic solution (plus signs) and the numerical solution after 30 Grad-Rubin 
iterations for the test case with periodic boundary conditions. The right panel shows -Ediv f° r the 
analytic solution (plus signs) and the numerical solution after 50 Grad-Rubin iterations for the test 
case with closed boundary conditions. The dashed lines are power-law fits (~ N"') to E^ lv for the 
numerical solutions. The power-law indices are 7 = —2.0 for the data in the left panel, and 7 = —1.4 
for the data in the right panel. The power-law index for -Ediv f° r the analytic solutions in both panels 
is 7 = —2.0. The fits are not shown. 



5.2. Test Case with Closed Side Boundaries 



1000 



For the second test we use the parameters k = tt/L, A = 0.9ir/L, and ao = 1. The grid 
sizes are the same as in Section 15.11 and we apply 50 Grad-Rubin iterations starting 
from a potential field. 

The right panel of Figure [2] shows Ai3 avg over 50 Grad-Rubin iterations for the 
three grid sizes used. For all three cases AB avg decreases exponentially for roughly 40 
iterations before becoming approximately constant. The rate of convergence does not 
appear to depend strongly on the grid size. 

Figure [7] shows the field lines of the analytic solution, the potential field, and the 
Grad-Rubin solution. The view in the two panels in the figure is along the y axis. The 
left panel shows the field lines of the initial potential field (blue field lines) and of the 
analytic solution (red field lines) . The right panel shows the field lines of the Grad-Rubin 
solution after 50 Grad-Rubin iterations (blue field lines) and of the analytic solution 
(red field lines). The Grad-Rubin solution closely matches the analytic solution. 

The right panel of Figure |4] shows the error metrics E m and Ecs for the initial 
potential field (plus signs and crosses), and for the numerical solution after 50 Grad- 
Rubin iterations (diamonds and squares). The error associated with the potential field 
is independent of N. For the Grad-Rubin solution both E m and Ecs decrease as power- 
laws, with indices 7 = —1.9 for E m and 7 = —2.0 for Ecs (based on least-squares fits, 
shown by the dashed lines). Both error metrics are found to scale approximately as 
I/TV 2 . The expected scaling, discussed in Section 15.2.21 is ~ 1/7V. The improvement 
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Figure 6. Execution time for 30 Grad- Rubin iterations as a function of grid size N for the first test 
case (Section 15,11 , The execution time has power-law scaling with a power-law index 7 = 3.8 based 
on a least-squares fit to the data (the dashed line). 

in the observed scaling may be due to the simplicity of the solution (for this solution 
By = 0) and is unlikely to occur more generally. 

The right panel of Figure [5] shows the error metric E^v as a function of N for 
the numerical solution after 50 Grad-Rubin iterations (diamonds) and for the analytic 
solution (plus signs). In this case i^div for the Grad-Rubin solution has power-law with 
index 7 = —1.0 (based on a least squares fit to the data, shown as a dashed line 
in the figure). The analytic solution has a power-law scaling in E^v with 7 = —2.0 
(the fit is not shown). There is a substantial difference between E^v for the analytic 
and numerical solution indicating that E^iv is providing a measure of the residual 
divergence of the numerical solution and is not due to the truncation error incurred by 
the numerical approximation to differentiation. 

6. Discussion and Conclusion 

We present an implementation of the Grad-Rubin method for solving the magnctostatic 
equations in a finite domain. The code allows two possible choices of boundary condi- 
tions on the side boundaries of the domain: either B is periodic on the side boundaries, 
or the normal component of B is zero on the side boundaries. We refer to these choices 
as periodic boundary conditions and closed boundary conditions. In both cases the 
normal component of B is zero on the top boundary of the domain. 

The code is tested in application to a simple analytic solution, for two choices of pa- 
rameters for the solution, which illustrate the periodic and closed boundary conditions 
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Figure 7. Comparison between the field lines for the analytic solution, the Grad-Rubin solution and 
the initial potential field for the test case in Section 15,21 The left panel shows the analytic solution 
(red field lines) and the potential field (blue field lines). The right panel shows the analytic solution 
(red field lines) and the numerical solution after 50 Grad-Rubin iterations (blue field lines). The view 
is along the y axis in both panels. The Grad-Rubin solution closely matches the analytic solution. 



respectively. In both cases the code accurately reconstructs the test solution. Several 
runs are performed for each test case with varying grid sizes, to demonstrate the scaling 
of the method with the size of the problem. The test case is highly idealized. The lower 
boundary conditions are exactly consistent with the magnetostatic model because they 
are derived from an exact solution, and the top and side boundary conditions exactly 
match the assumptions adopted for the numerical method. The idealization allows a 
rigorous test of the correctness of the implementation. 

For both test cases the Grad-Rubin method converges. The convergence is measured 
using the average absolute change AB avg in the field (see Equation ([iijl). The test 
case with periodic boundary conditions converges faster than the test case with closed 
boundary conditions, measured in the number of iterations. This is likely due to the 
second test case being significantly more non-potential than the first, rather than being 
due to the different boundary conditions. 

The code is found to accurately reproduce the analytic test cases. This is confirmed 
by a visual comparison of the field lines of the analytic solution and the reconstructed 
solution. For both choices of parameters/boundary conditions the numerical solution 
succeeds, based on this test. We also measure the success using the mean vector error 
E m which is defined by Equation (|45[) . and the Cauchy-Schwartz inequality based metric 
Ecs, which is defined by Equation (|4"6"]l . We also compute a measure of the residual 
divergence -Ediv which is defined by Equation (|47|) . In the absence of numerical error we 
would expect E nl = Eqs = E^w = 0. For both test cases E m , Ecs, and E^w decrease 
as N increases, with a power-law scaling in N, and we estimate the power-law index in 
each case from a least-squares fit to the data. The power-law indices are summarized in 
Table [TJ We obtain different scalings for the different metrics in each case. This may be 
attributed to the difference between the two test cases, and the different metrics. The 
slowest scaling achieved is ~ 1/N in each case, which is consistent with other Grad- 



Rubin implementations (jWheatlandl . 120061 ) . and with the estimate made in Section 
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13.2.21 The fastest scaling achieved is ~ 1/-/V 2 , which may be attributed to the simple 
form of the analytic solutions and is unlikely to be achieved for more general solutions. 
In general we expect a scaling ~ l/N. 

Table 1. Power law indices 7 for 
the scaling of the error metrics _E m 
, EcSi an d i?div with grid size. 



Test case £ m Eqs E div 

Periodic -1.3 -2.3 -2.0 
Closed -1.9 -2.0 -1.4 



The execution time of the code is found to sc ale as ~ N 4 , which is c omparable 
to the fastest force-free Grad-Rubin methods (e.g. IWheatlandl 20061 2007 ). However, 
the magnctostatic implementation is significantly slower in absolute terms than the 
force-free methods. The slowest step in both cases (force-free and magnetostatic) is the 
field line tracing used to update quantities in the volume at each iteration, and the 
magnctostatic case has two field line tracing steps (for and p( k >) compared to one 
(for a^) in the force-free case. 

This paper demonstrates the new method and code in application to a simple test 
case. The long-term goal of our work is to develop a method and code for reconstructing 
coronal magnetic fields for real solar data. Several obstacles remain to be overcome 
before this can be achieved. One problem is observational: the photospheric pressure 
profile is not currently available for real solar cases. However, progress might be made by 
assuming a simple model for the pressure, e.g. a constant /3 model, with p cx |B| 2 . Our 
model is also overly simplified in that we neglect gravity, and include only a pressure 
force. We are presently considering ways to include gravity also. Finally, it remains to be 
seen how well the techniques being developed cope with noisy data, and with data which 
are not strictly consistent with the magnetostatic model. Despite these qualifications, 
the method outlined here, and its successful application to analytic test cases, represents 
an important first step towards a Grad-Rubin method for magnetostatic reconstructions 
of the coronal magnetic field. 
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Appendix A 

A potential field Bo satisfying 

V x B = (48) 
and the boundary conditions at z = 

B z = Sobs (49) 

is used to initiate the Grad-Rubin iteration. The appropriate field with periodic side 
boundary conditions is given in Equation ([24"]) - ([26| in the text. For the choice of closed 
side boundaries the appropriate choice is 

oo oo 

B 0x (x,y,z) = c mn k m cosh(fc[z - L]) sin(fc m x) cos(k n y), (50) 

m— n—0 
oo oo 

B Qy (x,y,z) = c mn k n cosh(fc[z - L\) cos(fc m x) sin(fc»y), (51) 



m— n—0 



and 

oo oo 

B 0z (x,y,z) = -} ^2 c„ m fcsinh(A:[z - L}) cos(k m x) cos(fc„y), (52) 

m— n—0 

where k rn = 7rm/L, k n = nn/L, k 2 = fc^ + A;^, and where the coefficients c mn are given 

by 

Cmn^jzl [ dxdyB ohs (x, y)cos(k m x)coa(k n y)/ smh(kL). (53) 



o Jo 

3 „ A ;„tc ir, ~, AT3- 



Equations ([B1I|) - ([5^1) may be evaluated on a gr id with N 3 poin ts in ~ N 3 \og(N) 
operations using fast sine and cosine transforms ( Poularikas . 1996f ). 
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Appendix B 



This appendix presents the solution for the current-carrying component of the test field, 
B c , with closed top and side boundary conditions. The magnetic field may be calculated 
from a vector potential A in the Coulomb gauge (V • A = ) using V x A = B c . The 



vector potential satisfies the vector Poisson equation (jJacksonl 1199811 



V 2 A 



The boundary conditions for B c on the six plane boundaries are 

B c n = 0, 



(54) 



(55) 



where n is the unit vector normal to the boundary. Th e corresponding boundary condi- 
tions for the vector potential in the Coulomb gauge are (|Amari. Boulmezaoud. and Mikid .1 
1999h : 



and 



d n A = 0, 



A* =0, 



(56) 



(57) 



where A t denotes the component of A transverse to the boundary, and d n A. denotes 
the normal derivative. 

The vector potential A satisfying the boundary conditions Equations (|56D - ([5"T)) can 
be written as a Fourier series: 



oo oo oo 



A x = ^^^alll ip cos(k m x)sm(k n y)sm(k p z), (58) 

m— n— p— 

oo oo oo 



and 



m— n— p—0 



oo oo oo 



A z = ^2 ^2 a mn P sin(A; m a;) sm(k n y) cos(k p z), (60) 

m— n— p—0 

where k m = Ttm/L, k n = nn/L, and k p = Ttm/L. In the following we solve the 
Poisson equation to find an expression for a mnp . The process is demonstrated for the 
A x component, but the approach is similar for the other components. 

Substituting Equation ([55]) into the Poisson equation (Equation (|54"10 gives 



oo oo oo 



A x = 2J ^2 a m«P fc2 cos ( k mx) sm(k n y) sin(k p z) = (j,qJ x (x, y, z), (61) 

m— n— p—0 



where k 2 = fc^ + k 2 + k p . 

Applying the standard orthogonality relations (where 5 mn is the Kronecker delta): 

l L 
sin(7rms/L) sm(wns/L)ds = —S mn , (62) 

L L 

coa(iTms / L) cos(Trns/L)ds = —5 mn , (63) 
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and 



sm{iims/L) cos(irns/L)ds = 



(64) 



yields 



a 



(i) 

mnp 



8^0 
L 3 



L pL r L 



J x {x, y, z) cos(k m x) sin(fc„j/) sm(k p z)dxdydz. (65) 



'o Jo Jo 

Similar expressions apply for the other two coefficients: 

,-L P L pL 



a 



(2) 



and 



inn p 



8^o 
L 3 



JO JO 



J y {x, y, z) sin(k m x) cos(fc„y) sm(k p z)dxdydz, (66) 



8^o 
L 3 



L r L pL 



J z (x, y, z) sin(fc m a;) sin(fc„y) cos{k p z)dxdydz . (67) 

10 Jo Jo 

The magnetic field is obtained by evaluating B c = V x A. The components of B c are 



oo oo oo 

B C x{x, y, Z) = ^n a m!i P — k p 

m—0 n—0 p—0 

oo oo oo 

B C y{x, y, z) = [^p a mnp — ^ 

m—0 n—0 p—0 



(2) 
mnp 



m "'mnp 



sin(fc m x) cos(k n y) cos(k p z)/k 2 , (68) 

cos(k m x) sin(k n y) cos(/c p z)/fc 2 ,(69) 

(70) 



and 



oo oo oo 



B cz (x,y,z) = ^ X] X! [ fcma ™«P ~ fc " a mn P cos(/c m a;) cos(k n y) sm{k p z)/k 2 . (71) 



m—0 n—0 p—0 



The solution given by Equations (|6"8")) - (f7'l"j) may be computed on a grid with TV 3 points 
in ~ N 3 log(iV) operations using a combination of fast sine and cosine transforms. 
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